Assessing trends in population size of three unmarked species: A comparison of a multi‐species N‐mixture model and random encounter models

Abstract Estimation of changes in abundances and densities is essential for the research, management, and conservation of animal populations. Recently, technological advances have facilitated the surveillance of animal populations through the adoption of passive sensors, such as camera traps (CT). Several methods, including the random encounter model (REM), have been developed for estimating densities of unmarked populations but require additional information. Hierarchical abundance models, such as the N‐mixture model (NMM), can estimate abundances without performing additional fieldwork but do not explicitly estimate the area effectively sampled. This obscures the interpretation of its densities and requires its users to focus on relative measures of abundance instead. Hence, the main objective of our study is to evaluate if REM and NMM yield consistent results qualitatively. Therefore, we compare relative trends: (i) between species, (ii) between years and (iii) across years obtained from annual density/abundance estimates of three species (fox, wild boar and red deer) in central Spain monitored by a camera trapping network for five consecutive winter periods. We reveal that NMM and REM provided density estimates in the same order of magnitude for wild boar, but not for foxes and red deer. Assuming a Poisson detection process in the NMM was important to control for inflation of abundance estimates for frequently detected species. Both methods consistently ranked density/abundance across species (between species trend), but did not always agree on relative ranks of yearly estimates within a single population (between years trend), nor on its linear population trends across years (across years trend). Our results suggest that relative trends are generally consistent when the range of variability is large, but can become inconsistent when the range of variability is smaller.


| INTRODUC TI ON
Obtaining reliable demographic parameters, including (age-specific) survival, immigration, fecundity, and population size, is essential in wildlife management (Carpio et al., 2021;Williams et al., 2002).Since the estimation of population size over time enables population monitoring and is cheap in terms of data requirements (i.e., counts from point surveys replicated in space and time are sufficient), it is the (main) parameter of interest in many ecological studies.Moreover, collecting population counts has become increasingly cost-efficient over the last decades due to the adoption of automated sensor networks, such as camera traps (CT).Analytical frameworks for counts can yield precise estimates of population size (Keever et al., 2017;Palencia, Rowcliffe, et al., 2021) or trends therein (Kéry & Andrew Royle, 2010).When individuals are unmarked (i.e., they cannot be uniquely identified), obtaining population size using CTs, has been achieved through different analytical frameworks, including time-or space-to-event models (Moeller et al., 2018), distance sampling (Howe et al., 2017), random encounter (and staying time) models (Nakashima et al., 2018;Rowcliffe et al., 2008), spatial capture-recapture (Chandler & Royle, 2013) and site-structured abundance models (Kéry & Royle, 2016).
Broadly, these methods can be divided into two groups: (G1) those that estimate density from detection frequency and (G2) those that model animal counts as a function of an abundance and a detection parameter (i.e., detection probability or distance from an activity centre) that are jointly estimated (Loonam et al., 2021).
Importantly, the state variables of interest (density D, and abundance N) are slightly different across methods.In G1, density D represents the expected number of individuals N at any instant in time and within the collective set of camera viewsheds (i.e.areas in front of the CT in which individuals can be detected) with a total area A. In G2, N refers to local or site-abundance (henceforth "abundance"), the number of individuals available for detection during a specific survey duration and at a specific camera location.
In this study, we will evaluate whether relative trends in abundance/density are consistent across these two paradigms of treating unmarked population counts.Specifically, we will focus on the random encounter model (REM) and N-mixture model (NMM) (Royle, 2004) as representatives of G1 and G2, respectively.
We choose these models, as they are the most widely adopted methods for estimating unmarked population size using CT data (Gilbert et al., 2020).
The REM estimates animal density from trapping rate (aggregated count across the survey period), the average size of the detection zone of the camera, and average movement speed of the population under study (Rowcliffe et al., 2008).Obtaining speeds of movement can be done by tagging individuals with GPS collars, but increases the cost of the study and usually leads to underestimation of movement (Rowcliffe et al., 2012;Sennhenn-Reulen et al., 2017).However, a method to estimate the average speed of movement from CT pictures has recently been developed (Palencia, Fernández-López, et al., 2021;Rowcliffe et al., 2016).
REM makes the following assumptions: (i) detections are independent of each other, (ii) cameras are placed randomly relative to animal movement, (iii) individuals move independently of each other, and (iv) the populations under study are closed relative to the entire study area (i.e., no changes in overall population size within the survey period).The REM was found to be robust against violations in the independence of detections (Hayashi & Iijima, 2022), but not against non-random placement of cameras relative to animal movement (Cusack et al., 2015).Moreover, REM is sensitive to biased movement speeds, as well as to the method used to estimate the range of the camera viewsheds (Santini et al., 2022).
Nonetheless, Palencia, Rowcliffe, et al. (2021) obtained similar density estimates from REM compared to other methods representing G1.Finally, REM in its current form does not accommodate the modelling of spatial heterogeneity in density.
NMMs are hierarchical models that estimate the abundance at each camera (or site) based on counts from replicated surveys within the survey period rather than directly arriving at animal density for the collective set of camera viewsheds, as is done by REM.Consequently, the NMM requires that the study area is divided into discrete sites in order to infer abundance.The model assumes that (i) false-positive detections do not occur, (ii) detections are independent of each other, (iii) each individual has the same probability of being detected, and (iv) the local population size does not change throughout the survey period.As abundances are typically biased when some or all of these model assumptions are violated (Barker et al., 2018;Fogarty & Fleishman, 2021;Kéry & Royle, 2016;Link et al., 2018), solutions have been proposed that involve elegant ways to relax these assumptions (Dail & Madsen, 2011;Martin et al., 2011).Here, we formulate an NMM for open populations (Dail & Madsen, 2011), with a beta-Poisson detection process, building on ideas in Gomez et al. (2017) and Kéry and Royle (2016).Together, these adjustments accommodate changes in abundance between years and to some extent the occurrence of double counts (i.e., counting an individual twice during a survey).Furthermore, they allow the sharing of information on the detection process between commonly observed and rare species in a community (Gomez et al., 2017;Yamaura et al., 2011Yamaura et al., , 2012Yamaura et al., , 2016)).with fences impermeable to ungulates, such that the movement of wild boar and red deer in and out of the area should be limited.While these fences were not explicitly designed to be a movement barrier for fox, they may also hamper fox movement to some extent.

| Camera trapping network and data
Within the study area, a CT network was deployed each winter from 2017-18 to 2021-22.Each of these winters, 20 cameras were installed at the intersections of a lattice grid (with a camera spacing of ~2 km), which were fixed across years (Figure 1; Table 1).During winter 2019-20, eight additional cameras were placed.This resulted in camera locations, which are, to the best of our knowledge, random relative to the movements of the three target species, i.e., fox, red deer and wild boar.These species were selected as they were the only ones that generated sufficient records for analysis by REM/ NMM.As NMMs require data collected at discrete sites, we superimposed a hexagonal grid layer on the study area resulting in 336 grid cells of 0.2255 km 2 .This grid cell size trades off the possibility to capture fine-scale spatial variation with the possibility of detecting individuals at multiple trapping locations.The number of trapping days varied between years.During the winters of  and 2019-20, we used Bushnell Trophy Aggressor cameras, while Reconyx Hyperfire 2 and Browning Strike Force cameras were deployed during the winters of 2020-21 and 2021-22, respectively (Table 1).All cameras were mounted on trees ~50 cm above the ground, facing North and parallel to the ground.None of the cameras was baited to lure animals or was placed preferentially next to a trail.All cameras were set to be operative all day, to record a burst of consecutive photos (rapid fire) at each activation, and with the minimum triggering interval between activations.Timely check-ups were performed to determine battery levels and to verify that the cameras were still operating.Groups of consecutive photos were aggregated into sequences, which were manually annotated and used for the analysis of density/abundance by both REM and NMM (Bollen et al., 2023).

F I G U R E 1
Map of the study area in Quintos de Mora with hexagonal grids.The colour scale represents the proportion of forest in each grid (white: low proportion; dark brown: high proportion).Camera locations that had a camera deployed every winter are indicated by dots, those with deployments only in the winter of 2019-20 by triangles.The inset map shows the study area within Spain and Castilla-La Mancha (red).

| Statistical models
We analysed CT data for three different species using the REM and an extension of the NMM.While the REM was applied to each of the species specific data independently for each year, a single spatially and temporally explicit beta-Poisson NMM was fitted to the joint species data of all years in the study period.

| Random encounter model
We estimated animal density (individuals/km 2 ) for each target species and year in our study period separately using REM.Fitting REM requires three sources of information: (i) encounter rate y ∕ t , the rate at which individuals of a population come into contact with a CT, (ii) the radius r and angle of the camera viewsheds and (iii) animal movement speed v, which we obtained following Palencia, Rowcliffe, et al. (2021).First, we obtained the encounter rate of a target species by dividing the number of encounters y (i.e., total number of individuals from independent sequences of pictures of the species) by the total survey effort t (i.e., the sum of durations that each camera was active in the field).Note that we regarded pictures generated from each movement in and out of the camera viewshed as independent.Second, we estimated the effective radius r and angle ̂ by applying distance sampling to recorded positions (radial distance and angle) of each individual entering a camera viewshed (Rowcliffe et al., 2011).Third, we approximated the movement speed of each individual by dividing its distance travelled through the camera viewsheds by the time it took (i.e.time between first and last photo).For each target species, we then identified its main behavioural modes and averaged across all speed measurements of the corresponding mode to obtain behaviour specific speeds for the population.Day range was obtained by summation of the products of behaviour specific speeds and the proportion of time spent on each behaviour (Palencia, Fernández-López, et al., 2021).After obtaining encounter rate y ∕ t, the radius r and angle of the camera viewsheds, and the day range (movement speed) v for each population, animal densities were estimated by: Finally, we calculated standard errors associated with density estimates for each population using the delta method (Seber, 1982).

| N-mixture model
We modelled abundances for all target species and years in our study period using a single NMM.The NMM uses species counts, which are collected repeatedly in space and time, to model the expected number of individuals per site during a given survey period (abundance), knowing that some of the individuals that are present will not be detected (i.e., they do not show up in the species counts).
However, the discrete sites (0.2255 km 2 ) as we have defined them in section 2.2 are smaller than the typical home ranges of our target species.Hence, we cannot rule out that some individuals have been detected at multiple camera locations, violating the closure assumption.Thus, we interpreted abundances obtained from our NMM as relative abundances, i.e., the number of individuals that have used a site at least once during the survey period (henceforth "abundance" will refer to relative abundance).
We obtained species counts per survey day (24-h) by summation across all the individuals of that species counted on sequences of pictures from that day.This yielded counts y sijt for species s = 1, 2, 3 (fox, wild boar, red deer) at the subset of sites i = 1, 2, … , R that contained a CT during day j = 1, 2, … , J in year t = 1, 2, … , T. To correct for detection error NMM simultaneously estimates the detection probability (or rate) and the abundance of a species from y sijt .We fitted our NMM to species counts y sijt within a Bayesian estimation framework using Stan (via the R package cmdstanr), a probabilistic programming language that enables Bayesian estimation through a dynamic Hamiltonian Monte Carlo (HMC) sampler (Carpenter et al., 2017).Specifically, we assumed y sijt to be i.i.d.Poisson random variables, such that the detection process is given by: where the mean is a product of the latent number of individuals of species s at site i during year t (N sit ) and the species specific detection/trapping rate per camera per day (p sijt ).Note that by assuming a Poisson detection process, our NMM accommodates, to some extent, for double counts (i.e., it cannot account for individuals that are, on average, being detected >1 per survey day).We assumed that the trapping rate was constant across J days of year t and across R sites and that species specific detection rates p st followed a beta distribution: where we parameterised Beta( , ), such that = p t and = 1 − p t .
Under this parameterisation, p t and have a clear biological interpretation as the mean detection rate of all species in the community, and a measurement of similarity in species specific detection rates respectively (Dorazio et al., 2013).Furthermore, we modelled abundances N sit as a Poisson process with mean sit : Without further restrictions, the likelihood of this model involves an infinite sum over N sit , which we needed to restrict in order to sample from it.Therefore, we set species specific finite upper bounds which are much larger than the expected local population size ensuring that parameter estimates do not change appreciably beyond K s .
Moreover, we constructed the likelihood by marginalising over N sit 's with upper bound K given that Stan cannot sample discrete latent variables.
We defined two competing models,  1 and  2 , for which the detection process is identical.Both models estimate the detection rate of the community as a smooth curve f across years: However, the abundance process has an additional parameter in  2 compared to  1 , capturing a linear trend in abundance across years for each species, where s,0 and s,1 represent species specific intercepts and slopes for the trends across years, f ′ s models species specific smooth curves (trend noise), and f HSGP s is a spatial random effect.Both f and f ′ s use an exact Gaussian process (GP) (Golding & Purse, 2016;Williams & Rasmussen, 2006).For computational efficiency, we used the Hilbertspace reduced rank Gaussian process (HSGP) approach to model f HSGP s (Riutort-Mayol et al., 2020;Solin & Särkkä, 2020).As the inclusion of species specific random effects markedly increases the number of parameters, possibly resulting in models that are too difficult to fit, we also tested non-spatial versions of  1 and  2 omitting f HSGP s .
Prior specifications and goodness-of-fit diagnostics are detailed in Appendix A. We fitted all models using two parallel MCMC chains with 10,000 iterations, which included 5000 iterations that were discarded as burn-in iterations; this always resulted in satisfactory convergence (Table A1), following the guidelines by Vehtari et al. (2021).After fitting  1 and  2 (and their non-spatial versions), we performed a model selection by comparing their approximate leave-one-out expected log predictive densities (LOO-ELPD) (Vehtari et al., 2017).For a comparison of the results from a beta-binomial NMM and beta-Poisson NMM, we refer to Appendix B.

| Population trends
After model selection, we tested whether relative trends between species were consistent across the models by fitting a linear regression for yearly density (REM) versus abundance (NMM) estimates.We then compared temporal trends in density/abundance, obtained by REM and NMM in three ways.First, we computed the correlation between the ranks of relative trends between years in both methods using Spearman's rank correlation test.Next, we assessed the similarity of the trajectories of yearly growth rates, i.e., and also computed their Pearson correlations.Lastly, we compared slopes in linear trends across yearly densities/abundances.This is simply the estimated parameter ̂ s,1 of  2 for the NMM.However, to obtain this slope for REM, we needed to fit a linear regression line through estimates of yearly density post-hoc (for reference, we also did this for the NMM).Finally, to assess the precision of parameter estimates, we compared the coefficient of variation (CV) between yearly abundance and density.

| Trapping effort
Throughout the study period, we retain data from 4296 24-h periods

| Random encounter model
Mean annual densities estimated through REM lie between 0.41-0.73individuals/km 2 for fox, 5.34-7.14 for wild boar, and 25.06-46.63for red deer (Table 2).We do not observe a consistent increase or decrease in yearly densities for any of the target species.Relevant interannual variation is observed in the encounter rate and day range in all of the species (Table 2).Only seven fox encounters are recorded during 2020, hence we could not estimate the fox density for that year.

| N-mixture model
Model  1 has the best predictive performance according to LOO-ELPD, closely followed by  2 (Table 3).However, the standard error on the ΔLOO-ELPD between these models is substantial.
Furthermore, LOO-ELPD suggests that the spatial models are more consistent with the data than their non-spatial counterparts.
Both  1 and  2 fit the best to the fox data, followed by those of the red deer and finally the wild boar specific data (Table B2).For the remainder of the paper, all results from the NMM are generated based on the top-ranking model ( 1 ).The posterior mean detection rate of the community decreases until 2020-21 and shows a slight increase from 2020-21 to 2021-22 (Figure C1).For all the years in our analysis, there is a fair amount of posterior uncertainty, judging from the 50%, 80%, and 95% Bayesian credible intervals (BCI) for the community detection rate.The fox specific detection rate mostly resembles the community detection rate in terms of its mean trend and posterior uncertainty (Figure C1).The mean detection rates for both wild boar and red deer are distinct from the trend of the community.BCIs are narrow for all years and species.
Smooth temporal effects reveal year-to-year fluctuations in abundance of similar magnitude, but different trends between red deer and wild boar (Figure C2).The yearly variation in fox is larger than both of these species and also has a different trend.Smooth spatial effects display different magnitudes for all species (Figure C3).
The trend of spatial effects for foxes is not correlated with any of the spatial trends of other species (Figure C4).However, the spatial trends for red deer and wild boar are positively correlated.Together, relevant spatiotemporal variations in abundances are observed (Figures C5-C7).

| Population trends
Yearly densities (REM) and abundances (NMM) cannot be directly compared on their absolute scales (Figure 2c), yet they still contain important information on the consistency of species and/or year rankings across REM/NMM.The relationship between density (REM) and abundance (NMM) for the three species is captured well by a linear model (R-squared: 0.9141; Figure 2b).Only the ranks of yearly densities/abundances for the entire community are significantly correlated between  1 and REM (Table 4).However, both models produce similar trajectories in growth rates for fox, but not for wild boar and red deer (Figure 2c; Table 4).Interestingly,  1 and REM are in agreement about the direction of linear trends in density/ abundance estimates obtained post-hoc, except for fox (Figure 2d).

TA B L E 2
Values of the parameters of the estimated random encounter model (REM) for each population, where y/t is the encounter rate; v, the average distance travelled by an individual during a day (day range); r, the radius of detection; and Ɵ, the angle of detection.

TA B L E 3 Model comparison according
to the leave-one-out expected logpredictive densities (higher is better).
The 95% BCI of trend slopes from  2 ( s,1 ) overlap zero in all species.We did not attempt to compare the precisions of linear trends as they were obtained from values that are on substantially different scales.

| DISCUSS ION
In this study, we compared trends (i) between species, (ii) between years and (iii) across years obtained from empirical populations of three target species based on two models: REM and NMM.We have focussed on trends rather than absolute state variables (density/abundance) for two main reasons.First, state variables are slightly different across REM (density) and NMM (abundance), and also rely on different characterisations of space and time (Gilbert et al., 2020;Loonam et al., 2021).Second, the use of absolute population size in conservation and management has been challenged (Morellet et al., 2007), particularly when these are obtained through NMMs (Dennis et al., 2015;Gilbert et al., 2020).Since all of our results are based on empirical data, i.e., the truth is not known, we will focus our discussion on the consistency and precision of estimated trends rather than discussing their accuracy.Moreover, we note that in some cases our study may be underpowered to detect (small) differences in abundance given the modest number of cameras that we deploy, i.e., 20 (+8 in 2019-20).
Simulations may help to determine the number of cameras required to characterise trends in abundance (Ficetola et al., 2018).As the NMM is very sensitive to model assumptions, we tried to control for two common sources of bias in abundance: accidental double counting of unmarked species by a Poisson detection process (Kéry & Royle, 2016;Link et al., 2018), and unmodelled heterogeneity (Duarte et al., 2018;Link et al., 2018;Veech et al., 2016) by the inclusion of several fixed (Figure 2d) and random effects (Figures C2   and C3).Our study used different CT types across years, which has likely induced variability in the probability of detecting individuals that cross a camera viewshed.To account for this source of variability and potentially other differences leading to interannual variation in detectability (Hofmeester et al., 2019), we have included temporal effects in NMM and estimated angles and radii r separately for each year when applying REM.Moreover, we attempted to control bias that may result from inaccurate estimates of movement speed in REM, by correcting for different movement speeds in the main behavioural modes of a species (Palencia, Fernández-López, et al., 2021).
Although we warn users against using abundances from NMMs as absolute quantities (Barker et al., 2018), the yearly total abundances and densities retrieved from our NMM can be found in Table C1 for comparison with REM densities (Table 2).NMM and REM treat CT data differently (replicated counts vs. aggregated counts across the entire survey period).It is unclear how this impacts quantities derived from their state estimates, including relative trends, which is a study limitation.Regardless of the quantity of interest (either abundance or density), reaching a CV <0.25 is considered the minimum threshold for estimates to be useful for wildlife management (Skalski et al., 2005).The abundance estimates of our NMM meet this requirement across all species and years, while REM fails to deliver CV ≤ 0.25 for most species-year combinations (Table C2).However, some caution is warranted as higher precisions could result from overconfidence, rather than from correctly characterised improvements (Goldstein & De Valpine, 2022).
Possibly, the absence of movement parameters in the NMM, or the separation of model uncertainty over two subprocesses may lead to overconfidence in the precision of its abundances.On the contrary, underestimation of movement speed can lead to decreased precision of REM densities (Santini et al., 2022).Finally, weakly informative priors may have contributed to a lower CV (higher precision) in NMM abundances compared to densities from REM, which does not use priors.
Relative trends between species, based on species-rankings of density/abundance, over a 5-year period are consistent between NMM and REM (Figure 2b).A similar consistency in species-rankings based on relative abundance indices from camera surveys and densities from faeces counts has been observed (Ferretti et al., 2023).Species specific spearman rank correlations between yearly NMM abundances and REM reveal that these models capture different relative trends between years.This is reinforced by the differential progressions of the yearly growth rates for wild boar and red deer (Figure 2c).The reduced consistency in yearly growth rates relative to species-rankings between NMM and REM suggest that camera-based estimators will converge on similar qualitative results when the range of variability is large (i.e.between-species variability), but may not do so when the range of variability is smaller (i.e.between-year variability).However, fox growth rates reveal a similar progression between the models and are highly correlated showing that camera-based estimators may also yield consistent rank-order patterns when variability is smaller (Figure 2c; Table 4).Simulation studies should be performed to identify the exact circumstances resulting in consistent rank-order patterns between various camera-based estimators, and those resulting in a lack of consistency.Linear trends in density/abundance across years, obtained post hoc, were in agreement only for two of the three species (Figure 2d).However,  1 outperforming  2 in terms of LOO-ELPD suggests that these trends are stationary.
This is supported by the BCIs of s,1 overlapping zero.The precision of trend estimates was higher in NMM compared to REM for all species.This could be a consequence of our NMM, but not our REM, explicitly modelling a temporal dependency between yearly abundances.Modelling temporal autocorrelation can indeed shrink outlying observations to the mean trend (Outhwaite et al., 2018;Williams & Rasmussen, 2006).The possibility to model (non)linear trends in abundance, rather than obtaining them post hoc, is a main advantage of the NMM as it allows for hypothesis-testing through model selection (Table 3).This is also true for trends in space, which may help managers gain insight into the drivers of space use within a study area.In Quintos de Mora, for example, foxes seem to be more abundant in the west of the study area (Figure C5).trends, it holds the potential to be extended to a fully spatiotemporal analysis tool following ideas in Jousimo and Ovaskainen (2016).
Another possible avenue to capture spatial variation in density with REM, is to perform the analysis over a number of ecologically relevant strata.However, this may rapidly result in sample sizes that are too low to make meaningful inference when the number of strata increases or when the number of cameras available is limited.
In summary, we would advise practitioners against the use of NMMs when absolute densities/abundances are desired (although they may produce sensible estimates in some settings) and resort to REM or other methods instead (Efford & Fewster, 2013;Howe et al., 2017;Moeller et al., 2018;Nakashima et al., 2018).
However, practitioners should be aware that all of these methods require auxiliary data and/or field procedures, some of which can be time-consuming (Palencia, Rowcliffe, et al., 2021).Thus, when all that is needed are relative population trends, applying NMM may  B2), but yielded speedups in computation time (Table B1).Specifying a (beta-)Poisson instead of a (beta-)binomial model was effective against inflation of abundance for wild boar and red deer (Figure B2), and has probably helped to control bias from accidental double counting.

ADDITIONAL FIGURES AND TABLES FOR  1
In the main paper, we do not include visualisations for the detection rates, temporal and spatial random effects for the beta-Poisson NMM.
Here we show these results based on  1 , and we also include maps for the 2.5th and 97.5th percentiles of abundances.Furthermore, we give an overview of the posterior total abundance and density per species per year derived from  1 (Table C1).We computed the total abundance, and density based on  1 by: (i) calculating the expected number of individuals E N sit of species s at site i during year t, (ii) summing these expectations over all sites, such that E , and (iii) scaling them by the total surface area in the grid layer (i.e., the product of R , the number of sites, and A i , the area (in km 2 ) per site), hence E D st = E N st ∕ R ⋅ A i .Finally, we show comparisons of the coefficient of variation (CV) between abundances from  1 and densities from REM (Table C2).
Since the state variables of interest(density D and abundance N )    are different between REM and NMM, comparing these methods based on absolute estimates of their state variables would not yield T A X O N O M Y C L A S S I F I C A T I O N Population ecology, Spatial ecology valuable insights.Nevertheless, relative trends in their state variables should be largely consistent.However, rank-order estimates between REM and NMM may diverge in some circumstances given their different treatments of animal counts.Therefore, the objective of this study is to evaluate if REM and NMM yield consistent (relative) trends: (i) between species, (ii) between years and (iii) across years obtained from annual density/abundance estimates.Specifically, we evaluate i-iii empirically by fitting REM and NMM to CT data of fox (Vulpes vulpes), red deer (Cervus elaphus), and wild boar (Sus scrofa) from a Mediterranean area in central Spain collected during five consecutive winter periods.We believe that this comparison is relevant given the importance of population trends in wildlife conservation and management (see Prowse et al., 2021 for a recent example), and because it compares methods representing two fundamentally different paradigms of unmarked abundance estimation. 2 | MATERIAL S AND ME THODS 2.1 | Study area The study area (longitudes: 4.148-4.048°W; latitudes: 39.342-39.460°N) is the Quintos de Mora National Reserve.It has a total surface area of 68.64 km 2 , and is located south of the Montes de Toledo.The centre of the area is characterised by an open savanna, while the mountain ranges in the north and south are dominated by Mediterranean shrubland and natural forests (Figure 1).In the savanna, the most predominant species is Pinus pinea, while forests and shrubland are mainly composed of a mixture of Quercus coccifera, Quercus suber, Quercus ilex, Arbutus unedo, Erica spp, and Cistus spp.Quintos de Mora has altitudes ranging from 720 to 1050 m above sea level.The climate is slightly continental, characterised by cold winters and hot summers.Quintos de Mora has an annual precipitation between 300 and 400 mm.The entire study area is fenced

(
fox and wild boar) and 2189 24-h periods (red deer).This results in 2721, 520 and 226 observations of red deer, wild boar, and fox, respectively.The sampling period for fox and wild boar is extended relative to that of red deer, due to lower sample sizes in those species.Due to a defective camera, we retain data from only 19 CTs during the winters of2017-18, 2018-19, and 2020-21.

F
Consistency between population trends.(a) Mean ± 95% (B)CI abundances (NMM: closed circles ± full lines) and densities (REM: open circles ± dashed lines).(b) Linear trend ±95% CI bands for density (REM) versus abundance (NMM) estimates.95% (B)CI are displayed for each pair of estimates.(c) Mean ± 95% (B)CI growth rates in abundance (NMM: closed circles ± full lines) and in density (REM: open circles.95% CIs were not obtained due to a high error propagation using the delta method).(d) Mean ± 95% (B)CI slope coefficients for linear trends in yearly abundance (NMM) and density (REM) obtained by least squares regression.Linear trend in abundance captured by s,1 ( 2 ).No trend was visualised for  1 , as this model assumed that linear trends in abundances were absent.Colour scale -C: growth rate > 1 (red) or < 1 (blue), D: slope > 0 (red) or < 0 (blue).Fox density, and hence growth rate, in year 4 (2020-21) was not estimated by REM due to the low sample size.
Wild boar appear in three clusters, and red deer are clearly more abundant in the central savanna-like landscape.All of the species seem to avoid the mountainous, forested areas in the north and south of the study area to some extent.Regardless of the exact spatial ecology of the target species in Quintos de Mora, being able to capture meaningful information on animal space use patterns in general is of high ecological importance and should be considered when choosing a modelling framework.Although REM for camera trapping data currently lacks the possibility to model spatial and/or temporal be faster, cheaper and can provide insight in spatiotemporal dynamics of abundance.To bridge this gap we encourage the further A PPEN D I X B RESULTS FOR ALTERNATIVE OBSERVATION MODELS OF  1 In the main paper, we formulated  1 as a beta-Poisson NMM.Here we show the most important findings obtained by changing the observation process of  1 to either a Poisson, beta-binomial, or binomial distribution.Note that the beta-binomial NMM models species specific data jointly, just like the beta-Poisson NMM, while the Poisson NMM and binomial NMM need to be fitted to data of each species separately.In our study, a joint approach using a beta distribution did not qualitatively change the inference from NMM (Figures B1 and

F
Posterior mean (and 50, 80 and 95% Bayesian credible intervals) of yearly detection rates for each species, obtained from  1 .F I G U R E C 2 Posterior mean (and 50, 80 and 95% Bayesian credible intervals) of the smooth temporal variation (log-scale) in abundance of each species, obtained from  1 .F I G U R E C 3 Posterior mean of the smooth spatial variation (log-scale) in abundance of each species, obtained from  1 .For visual clarity, different scales are applied to each of the panels.F I G U R E C 4 Correlations between the species specific HSGP coefficients sm , obtained from  1 .F I G U R E C 6 2.5th Percentile of posterior yearly abundance for each species, obtained from  1 .Regions with more than 35% of forest cover are enclosed by the black lines.Heat scale: abundance -number of individuals that have used a grid cell at least once during the survey period.F I G U R E C 5 Posterior mean of yearly abundance for each species and each grid cell (0.2255 km 2 ), obtained from  1 .Regions with more than 35% forest cover are enclosed by black lines.Heat scale: abundance -number of individuals that have used a grid cell at least once during the survey period.F I G U R E C 7 97.5thPercentile of posterior yearly abundance for each species, obtained from  1 .Regions with more than 35% of forest cover are enclosed by the black lines.Heat scale: abundance -number of individuals that have used a grid cell at least once during the annual survey period.
Description of the yearly camera trapping survey.
, 1 − p t τ , TA B L E 1 Posterior means and 95% credible intervals for species specific overdispersion diagnostic Ĉ obtained from  1 and  2 .Posterior mean (and 95% Bayesian credible intervals) of yearly detection rate/probability for each species, obtained from  1 , a betabinomial version of  1 , and single-species versions of  1 with a Poisson or a binomial observation process.
TA B L E A 2 TA B L E B 1 Time (in seconds) needed to fit  1 , beta-binomial version  1 , and single-species versions of  1 with a Poisson or a binomial observation process for all species in the community.F I G U R E B 1 F I G U R E B 2 Posterior mean (95% Bayesian credible intervals) of yearly abundances for each species, obtained from  1 , a beta-binomial version of  1 , and single-species versions of  1 with a Poisson or a binomial observation process.
Note: Data represent posterior means (± posterior standard deviations).Estimated beta-Poisson N-mixture model parameter values for each population, where the abundance st is estimated directly, and the total abundance E N st and density E D st are derived from the number of sites R, and their respective surface area A i .Means ± standard deviations in abundance st for NMM, and in density abundance D st for REM.Estimated model parameter values for each population.
TA B L E C 2